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Abstract 

The problem of the uniqueness of solutions during the evaluation of primary 
energy spectra in the knee region using an extensive air shower (EAS) data set and 
the EAS inverse approach is investigated. It is shown that the unfolding of primary 
energy spectra in the knee region leads to mutually compensative pseudo solutions. 
These solutions may be the reason for the observed disagreements in the elementary 
energy spectra of cosmic rays in the 1-100 PeV energy range obtained from different 
experiments. 
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1 Introduction 



The Extensive air shower (EAS) inverse approach to a problem of the pri- 
mary energy spectra reconstruction in the region of 1 — 100 PeV energies has 
been an essential tool in the past decade [1,2,3,4,5,6,7]. Basically, it follows 
from the high accuracies of recent experiments [8,9,10,11,12] and the availabil- 
ity of the EAS simulation code [13], which was developed in the framework 
of contemporary interaction models in order to compute the kernel functions 
of a corresponding integral equation set [6,11]. At the same time, the energy 
spectra of primary (H, He and Fe) nuclei obtained from the KASCADE ex- 
periment [6] using the EAS inverse approach disagree with the same data from 
the ongoing GAMMA experiment [11,12], where parameterization of the EAS 
inverse problem is used. 

Below, a peculiarity of the EAS inverse problem is investigated, and one of 
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the possible reasons for the observed disagreements between the energy spectra 
in [6] and [11] is considered in the framework of the SIBYLL [14] interaction 
model. 



The paper is organized as follows: In Section 2 the EAS inverse approach 
and the definition of the problem of uniqueness is described. It is shown, 
that the abundance of primary nuclear species leads to pseudo solutions for 
unfolded primary energy spectra. The existence and significance of the pseudo 
solutions are shown in Section 4. The pseudo solutions for primary energy 
spectra were obtained on the basis of simulation of KASCADE [6] shower 
spectra. The EAS simulation model is presented in Section 3. In Section 5 
the peculiarities of the pseudo solutions are discussed in comparison with the 
methodical errors of the KASCADE data. 



2 Problem of uniqueness 



The EAS inverse problem is ill-posed by definition and the unfolding of 
the corresponding integral equations does not ensure the uniqueness of the 
solutions. The regularized unfolding on the basis of a priori information on 
expected solutions (smoothness, monotony and non-negativity) in some cases 
can redefine the inverse problem [15] and provide the appropriate solutions. 
However, the expected singularities (e.g. knees) in the primary energy spectra 
at 10 15 — 10 16 eV may erroneously be smoothed by regularization algorithms 
and vice versa, be imitated by the unavoidable oscillations [15] of the solu- 
tions. Furthermore, the EAS inverse problem implies evaluations of at least 
two or more unknown primary energy spectra from the integral equation set 
of Fredholm kind [6,11,12]. These peculiarities have not been studied in detail 
and the problem of the uniqueness of solutions can limit the number of eval- 
uated spectra. 

Let }a{E) be the energy spectrum of a primary nucleus A over the atmo- 
sphere, Wa{x\E) be the probability density function describing the transfor- 
mation of A and E parameters of the primary nucleus to a measurable vector 
x. Then the EAS inverse problem, i.e. the reconstruction of the energy spec- 
tra of Na primary nuclei on the basis of the detected spectra F(x) of EAS 
parameters, is defined by the integral equation 




a Na 



(1) 



A=A 1 



2 



Evidently, if f Al ,...A NA {E) are the solutions of eq. (1), the functions f A (E) + 
gA^E) should also be the solutions of (1), provided equation 



]T / g A (E)W A (x\E) dE = 0(±AY) 



(2) 



A 



is satisfied for the given measurement errors AY(x) and for at least one of the 
combinations of the primary nuclei 



The number of combinations (3) stems from a possibility of the existence of 
a set of functions g A (E) = gi jA (E), . . . , g ijA (E) for each of the primary nuclei 
(A), which can independently satisfy eq. (2). 

For example, suppose that N A = 3. Let us denote / g i: A k (E)W Ak (E) dE by 
h,A k and, for simplicity, set the right-hand side of eq. (2) to 0. Then, following 
expression (3), we find ric = 7 independent combinations of eq. (2): I\, Ak = 
for k = 1, 2 and 3, I 2 , Al + h,A 2 = 0, h, Al + h,A 3 = 0, I 3 ,a 2 + h,A 3 = and 
h,Ai + h,A 2 + h,A 3 = with different gi tAk {E) functions. The measurement 
errors ±AY on the right-hand side of these equations can both increase and 
decrease the domains of g it A k {E) functions. 

One may call the set of functions gA{E) the pseudo functions with the cor- 
responding pseudo solutions (spectra) fA(E)+gA(E). The oscillating qa(E) = 
g± iA (E) functions at j — 1 are responsible for the first N A equations / gi :A (E)W A (x.\E) dE = 
0(±AY), A = Ai, . . . A Na , due to the positive-definite probability density 
function W A (E). The pseudo solutions f A (E) + g 1A (E) can be avoided by 
using iterative unfolding algorithms [6,15]. 

Additional sources of the pseudo solutions originate from the mutually 
compensative effects at j > 2: 



inherent to eq. (2) for arbitrary groups of k and m ^ k primary nuclei. 
Since there are no limitations on the types of the pseudo functions (except 
for f A (E) + g A (E) > 0) that would follow from expression (4), and the num- 
ber of possible combinations (3) rapidly increases with the number of evalu- 
ated primary spectra (N A ), the problem of the uniqueness of solutions may 
be insoluble for iV^ > 3. Moreover, the pseudo functions have to restrict the 
efficiency of unfolding energy spectra for N A ~ 2 — 3, because the unification of 
Z — 1, ... ,28 primary nuclei spectra into 2 — 3 nuclear species (e.g. light and 
heavy) inevitably increases the uncertainties of the kernel functions W A (E) 




(3) 



£ / g Ak (E)W Ak (x\E)dE^ ]T / g Am (E)W A Jx\E)dE (4) 
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and thereby also increases the domains of the pseudo functions. 

Notice, that the pseudo solutions will always appear in the iterative unfold- 
ing algorithms if the initial iterative values are varied within large intervals. At 
the same time, it is practically impossible to derive the pseudo functions from 
the unfolding of equations (1,2) due to a strong ill-posedness of the inverse 
problem. However, for a given set of the measurement errors AY(x) and the 
known kernel functions Wa(x|£') for A = A±, . . . A^ A primary nuclei, eq. (2) 
can be regularized by parametrization of the pseudo functions ^(a, (3, . . . \ E). 
The unknown parameters (a, j3, . . .) can be derived from a numerical solution 
of parametric eq. (2), and thereby one may also evaluate the parametrized 
pseudo functions Qa{E). 

Below (Section 3), an EAS simulation model for computing the kernel func- 
tion Wa(E) and replicating the KASCADE [6] EAS spectral errors AF(x) is 
considered. 



3 EAS simulation model 

The primary energy spectra obtained in the KASCADE experiment were 
derived on the basis of the detected 2-dimensional EAS size spectra V(x) = 
Y (N e , N^) and an iterative unfolding algorithm [15] for Na = 5 primary nu- 
clei [6]. Evidently, whether these solutions are unique or not depends on the 
significance of the arbitrary pseudo functions \qa(E)\ from eq. (2). 

We suppose that the convolution of the shower spectra Wa(N € , \E) at 
the observation level and corresponding measurement errors a(N e ), 0"(A M ) 
[1] are described by 2-dimensional log-normal distributions with parameters 
£ e = In N e (A,E), ^ = In Nn(A,E), a e (A,E), a^A, E) and correlation coef- 
ficients p e ^(A,E) between the shower size (lnA" e ) and the muon truncated 
size (lnA^). We tested this hypothesis by the \ 2 goodness-of-fit test using 
the CORSIKA(NKG) EAS simulation code [13] for the SIBYLL2.1 [14] in- 
teraction model, 4 kinds of primary nuclei (A = p, He,0, Fe), 5 energies 
(E = 1, 3.16, 10, 31.6, 100 PeV) and simulation samples for each of E and A: 
5000, 3000, 2000, 1500, 1000 respectively in - 18° zenith angular interval. 
The values of corresponding x 2 (^ 5 Ej)/ n d.f., = -4, j = 1, . . .5) were 
distributed randomly in the interval 0.5 — 1.4 for the measurement ranges of 
the KASCADE experiment (A" e , min = 6.3 • 10 4 and A^ min = 4 • 10 3 ) and the 
bin sizes A In A^ e , A In = 0.075. 

Notice, that the combined 2-dimensional log-normal distributions with pa- 
rameters cr e ,i{A, E) at lnA e < £ e , a et2 {A,E) at lnA e > £ e , a^i(A,E) at 
lnA M < ^ and a^2(A,E) at InN^ > more precisely (x 2 / n d.f. < 1-2) de- 
scribe the shower spectra Wa(N £ , N^E) in the tail regions. 

We performed an additional test of the log-normal fit of the Wa spectra 
using multiple correlation analysis for the shower parameters simulated by the 
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log-normal Wa(N € , N^E) probability density functions and shower parame- 
ters obtained from the CORSIKA EAS simulations at power-law primary en- 
ergy spectra (7 = —1.5) and equivalent abundances of primary nuclei. The cor- 
responding correlation coefficients were equal to p(\nE\ In N e , In iV M ) = 0.97, 
pQnA\ lnA^lnAg = 0.71, p(\n A,\n N e ) = -0.14 ± 0.01, p(\nAAnN^) = 
0.18 ± 0.01, and were in close agreement for both methods of N e and 
generations. 

We replicated the KASCADE 2-dimensional EAS size spectrum Y(N e , A^) 
(and corresponding AY ) by picking out N e and randomly from the 2- 
dimensional shower spectra Wa(N £ , N^\E) after randomly picking A and E 
parameters of a primary particle from the power-law energy spectra 



with a rigidity-dependent knee E k = Z ■ 2000TV, the sharpness parameter 
e = 3 and normalization of the all-particle spectrum / J2a Ja{E) dE = 1. The 
relative abundance of nuclei was arbitrarily chosen to be 0.3,0.45,0.15 and 
0.1 for primary H,He,0 and Fe nuclei respectively, which approximately 
conforms with the expected abundance from balloon and satellite data [16]. 
The mediate values of the parameters of the probability density function 
Wa{N c , N m \E) were estimated by the corresponding log-parabolic splines. 

The total number of simulated EAS events was set to 7 ■ 10 5 in order to 
replicate the corresponding statistical errors AY(N e , N^) of the KASCADE 
data. 



4 Pseudo solutions 

On the basis of the obtained estimations of AY (N e , N^) (Section 3) for 
the KASCADE experiment, we examined the uniqueness of unfolding (1) by 
% 2 -the minimization: 



where G(N e>i , N^j) represents the left-hand side of eq. (2) for 2 kinds of the 
empirical pseudo functions 




(5) 




(6) 




(7) 
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gA (E) = a A {{\nE-(3 A f + r lA )(— ) , (8) 

while g A (E)+ f A (E) > 0, otherwise g A (E) = —f A (E). The unknown a A , j3 A , 7a 
and i] A parameters in expressions (7,8) were derived from \ 2 minimization (6). 
The numbers of bins were I = 60 and J = 45 with the bin size A In N e , A In ~ 
0.1. 

In fact, the minimization of x 2 (6) for different representations (7,8) of the 
pseudo functions g A (E) provides a solution of the corresponding parametric 
eq. (2) with a zero right-hand side. To avoid the trivial solutions g A (E) = 
and reveal the domains of the pseudo functions, the values of some of the 
parameters were arbitrarily fixed during the minimization of \ 2 (6). The mag- 
nitudes of the fixed parameters were empirically determined via optimization 
of conditions Xmm/ n d.f. — 1 an d \g A (E)\ ~ f A (E) for the pseudo spectra with 
the fixed parameters. 

The true primary energy spectra f A (E) for A = H, He,0, Fe nuclei (5) 
and the all-particle energy spectrum Y, }a{E) (lines) along with the corre- 
sponding distorted (pseudo) spectra f A (E) + g A (E) (symbols) are presented 
in Fig. 1 respectively. The parameters of the pseudo functions (7) derived for 
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Fig. 1. Primary energy spectra f A (E) and the all-particle spectrum J2 Ia(E) 

for A = H, He, O, Fe nuclei (lines) and the corresponding pseudo solutions 
f A (E) + g A (E) for the pseudo function (7) (symbols). 



Xmm/ n d.f. — 1-08 (^d/. = 717) are presented in Table 1. 

The effect of the pseudo functions (8) on the resulting primary energy 
spectra is shown in Fig. 2. Evaluations of the corresponding parameters are 
presented in Table 2 for Xmm/ n d.f. — 1-1. The variations of the cubic power 
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Table 1 

Parameters a a (TeV -1 ) and ja of the pseudo function (7) for different primary 
nuclei A and E m = 100 TeV. 

A a A ■ 10 4 1A 

p 1.10 ±0.06 2.71 ±0.04 

He -1.80 (fixed) 2.60 (fixed) 

O 0.97 ±0.05 2.65 ±0.04 

Fe -0.50 (fixed) 2.90 (fixed) 




indices in expression (8) in the range of 2 — 5 lead to different types of pseudo 
solutions as well. 

It is clear from Figs. 1,2, that the contribution of the pseudo functions 
Qa{E) can be comparable and even significantly larger than the values of the 
true spectra /a(-E')- Moreover, the pseudo solutions lose both the slopes and 
the intensities of the spectra. At the same time, the all-particle spectra slightly 
depend on the contribution of the pseudo functions. 

The same results (Tables 1,2) were obtained using both the combined 2- 
dimensional log- normal representation of the shower spectra Wa(N b , N^\E) 
(Section 3) and the 3-dimensional (InE, In N e , In JV^) parabolic interpolations 
of corresponding probability density functions obtained by the CORSIKA 
code. 
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Table 2 

Parameters a a (TeV -1 ), ja and r\ of the pseudo function (8) for different primary 
nuclei A and E m = 1000 TeV. 



A 


a A ■ 10 4 


Pa 


VA 


P 


-9.00 (fixed) 


7.76 ±0.01 


(fixed) 


He 


0.044 ± 0.02 


13.2 ± 1.08 


169 ± 98 





-0.80 (fixed) 


8.47 ±0.05 


0.94 ±0.16 


Fe 


0.010 ± 0.002 


11.4 ±0.14 


50 (fixed) 



Table 3 

Parameters aA (TeV -1 ) and ea (TeV) of the pseudo function (9) for different pri- 
mary nuclei A and Eh = 3000 TeV. 



A 


a A ■ 100 


£a/£h 


P 


-3.0 (fixed) 


1 (fixed) 


He 


3.05 ±0.07 


1.03 ±0.01 


O 


-0.84 ±0.06 


1.08 ± 0.03 


Fe 


0.15 ±0.02 


1.29 ±0.10 



Evidently, the range of relatively large measurement errors A7(x) expands 
the domain of the pseudo functions. Contributions of the mutually compen- 
sative effects (eqs. 2,4) of the pseudo functions to the domain of the pseudo so- 
lutions were tested using a 10 times larger EAS simulation sample (n = 7- 10 6 ) 
and the pseudo functions with evident singularity: 

g A (E) = a A eA l (—)\ (9) 



where 5 = — 1 at E < sa and 5 = —7 at E > ea- The singularity of the pseudo 
function (9) for A = H was fixed at Eh = 3000 TeV and the scale factor an = 
—0.03. The remaining parameters for primary nuclei A = He, O, Fe were es- 
timated by x 2 -minimization (6) and presented in Table 3 for Xmm/ n d.f. = 2-01 
and n d _f, = 857. The accuracies of integrations (2) were about 0.1%. The cor- 
responding pseudo solutions are shown in Fig. 3. 

Since the measurement errors are negligibly small, the significance of the 
mutually compensative effects is well seen. The singularity of the proton spec- 
trum was approximately compensated by the He and O spectra. This is due to 
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Fig. 3. The same as Fig. 1 for the pseudo function (9) and n = 7 • 10 6 simulated 
showers. 



both the large number (jiq = 15) of possible mutually compensative combina- 
tions (3) and the peculiarities of EAS development in the atmosphere (kernel 
functions Wa(E), Section 3), which are expressed by the approximately log- 
linear dependences of the statistical parameters < lniV e >, < lniV M >, a e 
and <r M of shower spectra Wa{E) on energy (InE) and nucleon number (In A) 
of primary nuclei [20,21]. The value of xLm/ n d.f. f° r a 10 times smaller EAS 
sample (n — 7 • 10 5 ) was equal to 0.25. 



5 Discussion 

The results from Figs. 1-3 show that the pseudo functions with mutually 
compensative effects exist and belong practically to all families - linear (7), 
non-liner (8) and even singular (9) in a logarithmic scale. 

The all-particle energy spectra in Figs. 1-3 are practically indifferent to the 
pseudo solutions of elemental spectra. This fact directly follows from eq. (2) for 
pseudo solutions and is well confirmed by the identity of the GAMMA [11,12] 
and KASCADE [6] all-particle energy spectra in spite of disagreements of the 
elemental (p,He,Fe) primary energy spectra (see [11,12]). 

The x 2 minimization (6) uses mainly the nearest pseudo energy spectra with 
free parameters for compensation of the pseudo spectra with fixed parameters. 
The significance of the pseudo functions |g^(i?)| in most cases exceeds the sig- 
nificance of the evaluated primary energy spectra }a{E) and unfolding of (1) 
can not be effective for Na = 4. 
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Fig. 4. Domains of the pseudo solutions for He and Fe primary nuclei (light shaded 
areas) and corresponding "methodical errors" of the KASCADE unfolding spectra 
(dark shaded areas) taken from [6] . The solid and dotted lines resulted from pseudo 
functions (7) and the dashed lines stemmed from (8). 



The unfolding of the primary energy spectra for Na = 5 will increase the 
number of possible combinations (3) of the pseudo solutions and the corre- 
sponding pseudo functions by a factor of two. Taking into account the large 
values of applied Xmin/ ' n d.f. — 2 — 3 [6] one may conclude that the contribu- 
tions of the pseudo functions in the unfolded energy spectra of [6] have to be 
dominant. 

The "methodical errors" obtained in [6] for Na = 5 define the uncertainties 
of the solutions intrinsic only to the given unfolding algorithms. The existence 
and significance of the mutually compensative pseudo solutions follow from 
eqs. (1,2) and from the peculiarities of the shower spectra Wa(x\E) regardless 
of the unfolding algorithms. 

Comparison of the methodical errors (J 'a(E) + AJ 'a(E)) / J 'a(E) for A = He 
and A = Fe from [6] with corresponding errors (/a(-E') + Qa{E))/ Ja{E) due 
to the pseudo solutions from expressions (7,8) are shown in Fig. 4. The mag- 
nitudes of the fixed parameters were empirically determined by maximizing 
\gHe(E)\ (left panel) and \gpe{E)\ (right panel) for a given goodness-of-fit test 
X 2 mi n/ n d.f - 2.5 from [6]. 

It is seen that the methodical errors (dark shaded areas) from [6] signif- 
icantly underestimate the contribution of the pseudo solutions (light shaded 
areas) from expressions (7, 8). Moreover, the methodical errors from [6] slightly 
depend on the primary energy (or statistical errors), whereas the domains of 
the pseudo solutions strongly correlate with the statistical errors according to 
definition (2). 
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6 Conclusion 



The results show that the reconstruction of primary energy spectra using 
unfolding algorithms [6,15] can not be effective and the disagreement between 
the KASCADE [6] and GAMMA [11,12] data is insignificant in comparison 
with the large domains of the mutually compensative pseudo solutions (Fig 4) 
of the unfolded spectra [6]. 

Even though the oscillating pseudo solutions g^A(E) (Section 2) are possi- 
ble to avoid using regularization algorithms [15], the mutually compensative 
effect (4) of the arbitrary pseudo functions Qa{E) intrinsic to the expression 
(2) is practically impossible to avoid at Na > 1. 

The uncertainties of solutions due to the mutually compensative pseudo 
functions can be obtained by varying the initial values of iterations within a 
wide range in the frameworks of a given unfolding algorithm. 

To decrease the contributions of the mutually compensative pseudo solu- 
tions one may apply a parameterization of the integral equations (1) [1,2,4,11,12] 
using a priori (expected from theories [17,18,19]) known primary energy spec- 
tra with a set of free spectral parameters. This transforms the EAS inverse 
problem into a set of equations with unknown spectral parameters, and thereby 
the EAS inverse problem is transmuted into a test of the given primary energy 
spectra using detected EAS data [4]. The reliability of the solutions can be 
determined by their stability depending on the number of spectral parame- 
ters, the agreement between the expected and detected EAS data sets, and 
the conformity of the spectral parameters with theoretic predictions. 

The all-particle energy spectra (Fig. 1-3) are practically indifferent toward 
the pseudo solutions for elemental spectra. 

The obtained results depend slightly on the spectral representations of the 
shower spectra Wa{E) and the primary energy spectra /a(-E). 
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